function [ u ] = evaluateInteriorPoint( x, y, domain2D, theta )
%EVALUATEINTERIORPOINT Summary of this function goes here
%   Detailed explanation goes here

global magicE;

n = size(x, 1);
% nV = size(domain2D.vx, 1);
nE = size(domain2D.e1, 1);

u = zeros(n, 2);

epsilon = magicE;

for i = 1 : n
    
    p = [x(i), y(i)];
    
    for j = 1 : nE
      qi0 = domain2D.e1(j);
      qi1 = domain2D.e2(j);
      
      vqn0 = [domain2D.vnx(qi0), domain2D.vny(qi0)];
      vqn1 = [domain2D.vnx(qi1), domain2D.vny(qi1)];
      
      q0 = [domain2D.vx(qi0), domain2D.vy(qi0)] + epsilon * vqn0;
      q1 = [domain2D.vx(qi1), domain2D.vy(qi1)] + epsilon * vqn1; 
      
%       qn = [domain2D.nx(j), domain2D.ny(j)];
      
      v0 = -1 * evaluateQuadraturePoint(p, q0);
      v1 = -1 * evaluateQuadraturePoint(p, q1);
      
      ql = domain2D.l(j);
      
      u(i, :) = u(i, :) + ((v0 + v1) * ql / 2.0 * theta(2*j-1:2*j))'; 
      
%       A(2*i-1:2*i, 2*j-1:2*j) = A(2*i-1:2*i, 2*j-1:2*j) + (v0 + v1) * ql / 2;
      
   end
    
%     for j = 1 : m        
%         
%         qi0 = j;
%         qi1 = mod(j, m) + 1;
%         
%         q0 = [domain2D.x(qi0), domain2D.y(qi0)];
%         q1 = [domain2D.x(qi1), domain2D.y(qi1)];
%               
%         qn = [domain2D.nx(j), domain2D.ny(j)];
%         
%         v0 = -1 * evaluateQuadraturePoint(p, q0);
%         v1 = -1 * evaluateQuadraturePoint(p, q1);
%         
%         u(i, :) = u(i, :) + ((v0 + v1) * domain2D.l(j) / 2.0 * theta(2*j-1:2*j))'; 
%     end
    
end

end

